Introducción

En R podemos encontrar una serie de paquetes que se pueden usar para leer, visualizar y analizar datos espaciales. En esta ocasión nos centraremos en la manipulación de datos raster, específicamente nos centraremos en los paquetes:

if (!require(raster)) install.packages('raster',dependencies=T)
library(raster)
if (!require(rasterVis)) install.packages('rasterVis',dependencies=T)
library(rasterVis)
if (!require(velox)) install.packages('velox',dependencies=T)
library(velox)
if (!require(fasterize)) install.packages('fasterize',dependencies=T)
library(fasterize)
if (!require(sf)) install.packages('sf',dependencies=T)
library(sf)
if (!require(rgdal)) install.packages('rgdal',dependencies=T)
library(rgdal)
if (!require(tidyverse)) install.packages('tidyverse',dependencies=T)
library(tidyverse)

Iniciaremos cargando un modelo digital de elevación.

MDE<-raster("Datos/Michoacan_60m.tif")

Al impimimir el objecto RasterLayer se muestra el resumen de los datos

MDE
## class       : RasterLayer 
## dimensions  : 4463, 6615, 29522745  (nrow, ncol, ncell)
## resolution  : 0.000555556, 0.000555556  (x, y)
## extent      : -103.7381, -100.0631, 17.91491, 20.39435  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/alequech/Documents/GitHub/Charla_RSpatialES_raster-/Datos/Michoacan_60m.tif 
## names       : Michoacan_60m 
## values      : -32768, 32767  (min, max)

ahora realizemos el plot con de nuestro MDE

plot(MDE)

Corte y confección

El siguiente paso consiste en recortar el MDE para un area de interes especifica para esto emplearemos la funcion crop(). Para esto cargaremos primero un archivo vectorial usnado shapefile().

areaI<-shapefile("Datos/mun.shp")
areaI
## class       : SpatialPolygonsDataFrame 
## features    : 4 
## extent      : -102.5759, -101.9362, 19.17302, 19.62749  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       :               NOM_MUN, Id 
## min values  : Nuevo Parangaricutiro,  1 
## max values  :               Uruapan,  4
plot(areaI)

MDE_I<-crop(MDE,areaI)
MDE_I
## class       : RasterLayer 
## dimensions  : 818, 1152, 942336  (nrow, ncol, ncell)
## resolution  : 0.000555556, 0.000555556  (x, y)
## extent      : -102.5759, -101.9359, 19.17324, 19.62769  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /private/var/folders/z6/7nkbhk316q1d2s06j6kybk5h0000gn/T/RtmphqdBoy/raster/r_tmp_2018-11-29_205627_750_14741.grd 
## names       : Michoacan_60m 
## values      : 315, 3843  (min, max)
plot(MDE_I)

plot(MDE_I)
plot(areaI, add=T)

MDE_II<-mask(MDE_I,areaI)
MDE_II
## class       : RasterLayer 
## dimensions  : 818, 1152, 942336  (nrow, ncol, ncell)
## resolution  : 0.000555556, 0.000555556  (x, y)
## extent      : -102.5759, -101.9359, 19.17324, 19.62769  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : Michoacan_60m 
## values      : 600, 3843  (min, max)
plot(MDE_II)
plot(areaI,add=T)

Es posible explorar nuestro conjunto de datos usando hist y boxplot

hist(MDE_II)

boxplot(MDE_II)
## Warning in .local(x, ...): taking a sample of 1e+05 cells

Una operación bastante común cuando se trabaja con información espacial es el cambio de sistemas de referencia. Es decir proyectar de un sistema a otro para esto es crucial conocer el sistema de referencia de destino. Con el paquete raster podemos usarprojectRaster() para este fin.

En este caso voy a emplear el sistema de referencia Mexico ITRF2008 / LCC, para mayor informacion puede consultar http://geotiff.maptools.org/proj_list/

crs(MDE_II)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
?projectRaster()
Rcrs<-CRS("+init=epsg:4326")
Rcrs
## CRS arguments:
##  +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## +towgs84=0,0,0
#usando PROJ.4
RPROJ.4<-CRS("+proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") 
RPROJ.4
## CRS arguments:
##  +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102
## +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs
#MDE_LCC<-projectRaster(MDE_II,crs=Rcrs)
MDE_LCC<-projectRaster(MDE_II,crs=RPROJ.4)
MDE_LCC
## class       : RasterLayer 
## dimensions  : 830, 1164, 966120  (nrow, ncol, ncell)
## resolution  : 58.2, 61.3  (x, y)
## extent      : 2439291, 2507036, 797136.1, 848015.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : Michoacan_60m 
## values      : 600, 3833.119  (min, max)
pendiente<-terrain(MDE_LCC, opt='slope',unit='degrees')
plot(pendiente)

Es posible usar operadores lógicos para seleccionar un conjunto de pixeles que cumpla con un criterio determinado. Por ejemplo pixeles con un valor mayor a 10º.

consulta1 <- pendiente > 10
consulta1
## class       : RasterLayer 
## dimensions  : 830, 1164, 966120  (nrow, ncol, ncell)
## resolution  : 58.2, 61.3  (x, y)
## extent      : 2439291, 2507036, 797136.1, 848015.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : 0, 1  (min, max)
plot(consulta1, col = c("gray", "green"))

consulta2 <- pendiente 
consulta2[ pendiente < 10]<-NA
consulta2
## class       : RasterLayer 
## dimensions  : 830, 1164, 966120  (nrow, ncol, ncell)
## resolution  : 58.2, 61.3  (x, y)
## extent      : 2439291, 2507036, 797136.1, 848015.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : slope 
## values      : 10.00002, 53.83355  (min, max)
plot(consulta2)

consulta3 <- pendiente 
consulta3[ pendiente < 10 | pendiente > 40]<-NA
consulta3
## class       : RasterLayer 
## dimensions  : 830, 1164, 966120  (nrow, ncol, ncell)
## resolution  : 58.2, 61.3  (x, y)
## extent      : 2439291, 2507036, 797136.1, 848015.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : slope 
## values      : 10.00002, 39.99954  (min, max)
plot(consulta3)

El proceso de reclasificación es tarea cotidiana en el análisis espacial, en el siguiente ejemplo empleamos projectRaster(). En este caso

m <- c(0, 15, 1,  15, 20, 2,  20, 60, 3)
clases <- matrix(m, ncol=3, byrow=TRUE)
clases
##      [,1] [,2] [,3]
## [1,]    0   15    1
## [2,]   15   20    2
## [3,]   20   60    3
rec = reclassify(pendiente, clases)
rec
## class       : RasterLayer 
## dimensions  : 830, 1164, 966120  (nrow, ncol, ncell)
## resolution  : 58.2, 61.3  (x, y)
## extent      : 2439291, 2507036, 797136.1, 848015.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : 0, 3  (min, max)
plot(rec)

Agregación de datos espaciales

Es posible modificar el arreglo espacial de las celdas agrupando celdas. La función aggregate() permite reagrupar varios pixeles, obteniendo un mapa raster de menor resolución espacial. La función fun() permite controlar el cálculo del valor del pixel de baja resolución que corresponde a varias celdas en los datos de entrada. En el siguiente ejemplos se reagrupan 2x2 pixeles, se promedian los valores de las cuatro celdas.

MDE_agreg <- aggregate(MDE_LCC, fact=2, fun="mean")
MDE_agreg
## class       : RasterLayer 
## dimensions  : 415, 582, 241530  (nrow, ncol, ncell)
## resolution  : 116.4, 122.6  (x, y)
## extent      : 2439291, 2507036, 797136.1, 848015.1  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : Michoacan_60m 
## values      : 603.2935, 3818.155  (min, max)

Operaciones zonales

El paquete raster permite también llevar a cabo operaciones zonales que consiste en calcular algún indice sobre una capa estratificando el cálculo en función en otra capa.
Estadisticas zonales, Fuente: http://eujournal.org/files/journals/1/books/JeanFrancoisMas.pdf

areaI_LCC<- spTransform(areaI, RPROJ.4)
areaI_LCC.sp<- st_as_sf(areaI_LCC)
system.time(rasterize(areaI_LCC,MDE_LCC,field='Id'))
##    user  system elapsed 
##  26.486   0.151  26.761
system.time(fasterize(areaI_LCC.sp,MDE_LCC,field='Id'))
##    user  system elapsed 
##   0.011   0.005   0.030
mun_R<-fasterize(areaI_LCC.sp,MDE_LCC,field='Id')
plot(mun_R)

Elev_P= zonal(MDE_LCC, mun_R, fun = "mean") %>%
  as.data.frame()
Elev_P
##   zone     mean
## 1    1 2021.665
## 2    2 1618.826
## 3    3 2235.664
## 4    4 1906.949

Operaciones Focales

raster permite realizar operaciones de filtrado espacial. La operacion focal a continuación calcula desviacion estandar de los valores de las celdas en una ventana móvil de 7 x 7 pixeles.

matrix(1,7,7)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    1    1    1    1    1    1    1
## [2,]    1    1    1    1    1    1    1
## [3,]    1    1    1    1    1    1    1
## [4,]    1    1    1    1    1    1    1
## [5,]    1    1    1    1    1    1    1
## [6,]    1    1    1    1    1    1    1
## [7,]    1    1    1    1    1    1    1
#pendiente.sd <- focal(pendiente, w = matrix(1,7,7), fun = sd,pad=T)
pendiente.sd <- focal(pendiente, w = matrix(1,7,7), fun = sd)
plot(pendiente.sd)

Algebra de mapas

Para los siguientes ejemplos se usaran datos del producto MOD13Q1 de MODIS, específicamente la capa de NDVI. En este producto los valores de NDVI tienen un factor de escala de 0.0001.

modis1<-raster("Datos/Modis/MOD13Q1_NDVI_2010_001.tif")
modis1
## class       : RasterLayer 
## dimensions  : 2359, 1836, 4331124  (nrow, ncol, ncell)
## resolution  : 232, 234  (x, y)
## extent      : 43895.63, 469847.6, 2449977, 3001983  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/alequech/Documents/GitHub/Charla_RSpatialES_raster-/Datos/Modis/MOD13Q1_NDVI_2010_001.tif 
## names       : MOD13Q1_NDVI_2010_001 
## values      : -2000, 9487  (min, max)

Ahora procederemos a escalar los valores de NDVI

system.time(modis1/10000)
##    user  system elapsed 
##   0.352   0.053   0.471
system.time(calc(modis1,function(x){x /10000}))
##    user  system elapsed 
##   0.478   0.020   0.498

Los datos raster puede eventualmente manejarse como imágenes multibanda con stack() o brick().

lmodis<-list.files(pattern = "MOD13*", recursive = TRUE, include.dirs = TRUE)
lmodis
## [1] "Datos/Modis/MOD13Q1_NDVI_2010_001.tif"
## [2] "Datos/Modis/MOD13Q1_NDVI_2010_017.tif"
## [3] "Datos/Modis/MOD13Q1_NDVI_2010_033.tif"
## [4] "Datos/Modis/MOD13Q1_NDVI_2010_049.tif"
## [5] "Datos/Modis/MOD13Q1_NDVI_2010_065.tif"
## [6] "Datos/Modis/MOD13Q1_NDVI_2010_081.tif"
## [7] "Datos/Modis/MOD13Q1_NDVI_2010_097.tif"
modis2010<-stack(lmodis)
modis2010
## class       : RasterStack 
## dimensions  : 2359, 1836, 4331124, 7  (nrow, ncol, ncell, nlayers)
## resolution  : 232, 234  (x, y)
## extent      : 43895.63, 469847.6, 2449977, 3001983  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : MOD13Q1_NDVI_2010_001, MOD13Q1_NDVI_2010_017, MOD13Q1_NDVI_2010_033, MOD13Q1_NDVI_2010_049, MOD13Q1_NDVI_2010_065, MOD13Q1_NDVI_2010_081, MOD13Q1_NDVI_2010_097 
## min values  :                 -2000,                 -2000,                 -2000,                 -2000,                 -2000,                 -2000,                 -2000 
## max values  :                  9487,                  9818,                  9446,                  9812,                  9402,                  9216,                  9774
#modis2010R<-mosdis2010/100000
modis2010R<-calc(modis2010,fun=function(x){x / 10000})
modis2010R
## class       : RasterBrick 
## dimensions  : 2359, 1836, 4331124, 7  (nrow, ncol, ncell, nlayers)
## resolution  : 232, 234  (x, y)
## extent      : 43895.63, 469847.6, 2449977, 3001983  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /private/var/folders/z6/7nkbhk316q1d2s06j6kybk5h0000gn/T/RtmphqdBoy/raster/r_tmp_2018-11-29_205817_750_82109.grd 
## names       : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7 
## min values  :    -0.2,    -0.2,    -0.2,    -0.2,    -0.2,    -0.2,    -0.2 
## max values  :  0.9487,  0.9818,  0.9446,  0.9812,  0.9402,  0.9216,  0.9774

Como el usuario debe suponer en este momento también es posible realizar operaciones entre capas raster, nuevamente nos encontramos con diferentes maneras de abordar el problema.

r1<-raster(modis2010R,1)
r1
## class       : RasterLayer 
## band        : 1  (of  7  bands)
## dimensions  : 2359, 1836, 4331124  (nrow, ncol, ncell)
## resolution  : 232, 234  (x, y)
## extent      : 43895.63, 469847.6, 2449977, 3001983  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /private/var/folders/z6/7nkbhk316q1d2s06j6kybk5h0000gn/T/RtmphqdBoy/raster/r_tmp_2018-11-29_205817_750_82109.grd 
## names       : layer.1 
## values      : -0.2, 0.9487  (min, max)
r2<-modis2010R[[2]]
r2
## class       : RasterLayer 
## band        : 2  (of  7  bands)
## dimensions  : 2359, 1836, 4331124  (nrow, ncol, ncell)
## resolution  : 232, 234  (x, y)
## extent      : 43895.63, 469847.6, 2449977, 3001983  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /private/var/folders/z6/7nkbhk316q1d2s06j6kybk5h0000gn/T/RtmphqdBoy/raster/r_tmp_2018-11-29_205817_750_82109.grd 
## names       : layer.2 
## values      : -0.2, 0.9818  (min, max)
promedio<- (r1+r2) / 2
plot(promedio)

Por supuesto en este caso tenemos un conjunto de datos de 13 capas, en este caso podemos emplear la función overlay(), esta función es similar al paquete calc(). La diferencia radica en que overlay() esta pensado para agilizar las operaciones entre capas.

#NDVI_promedio<-mean(modis2010R)
NDVI_promedio<-overlay(modis2010R, fun=mean)
NDVI_promedio
## class       : RasterLayer 
## dimensions  : 2359, 1836, 4331124  (nrow, ncol, ncell)
## resolution  : 232, 234  (x, y)
## extent      : 43895.63, 469847.6, 2449977, 3001983  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /private/var/folders/z6/7nkbhk316q1d2s06j6kybk5h0000gn/T/RtmphqdBoy/raster/r_tmp_2018-11-29_205905_750_42759.grd 
## names       : layer 
## values      : -0.1297, 0.9007429  (min, max)
plot(NDVI_promedio, main="Promedio de NDVI de los primeros 100 dias del año 2010")

Para guardar un archivo raster es necesario conocer primero el tipo de dato que deseamos guardar. Para mayor información consultar (https://www.rdocumentation.org/packages/raster/versions/2.8-4/topics/dataType)

dataType(NDVI_promedio)
## [1] "FLT4S"
writeRaster(NDVI_promedio,'NDVI_P.tif', format="GTiff",datatype="FLT4S", overwrite=TRUE)

Mejorar la visualización

El paquete rasterVis define métodos adicionales para el acceso, la manipulación y la visualización de datos raster.

plot(MDE_LCC,
     main="Modelo Digital de Elevación",
     cex.main = 0.7,
     col=viridisLite::viridis(10)) 
plot(areaI_LCC, add=T)

Aquí aplicamos la función leveplot () que está optimizada para mostrar datos ráster y que proporciona algunas características gráficas adicionales.

m <- 100
levelplot(MDE_LCC,
          margin=FALSE,
          contour=TRUE,
          at=seq(0, maxValue(MDE_LCC), m),
          main = paste('Modelo Digital de Elevación, contorno', m, 'm'),
          par.settings=rasterTheme())

histogram(MDE_LCC)

levelplot(modis2010R)

Recomendaciones